*PATCH,GEXAM1.
*CMZ :  3.21/02 29/03/94  15.41.35  by  S.Giani
*DECK, GTCKOV
*CMZ :          12/09/95  11.07.14  by  S.Ravndal
*CMZ :  3.21/04 13/12/94  15.17.13  by  S.Giani
*-- Author :
      SUBROUTINE G3TCKOV
C.
C.    ******************************************************************
C.    *                                                                *
C.    *   This routine is called to follow the Cherenkov photons       *
C.    *   created during the tracking of charged particles and         *
C.    *   simulate the relevant processes along the way, until either  *
C.    *   the photon is absorbed or exits the detector. Processes      *
C.    *   currently simulated are absorption in-flight, and reflection *
C.    *   /transmission/absorption at a medium boundary. There are two *
C.    *   boundary types: dielectric-metal and dielectric-dielectric.  *
C.    *   For each of these there is a continuum of reflectivity       *
C.    *   and of surface quality from mirror finish to matte. The      *
C.    *   surface model is contained in routine GHSURF.                *
C.    *                                                                *
C.    *   ==>Called by : G3TRACK                                       *
C.    *      Authors    F.Carminati, R.Jones ************              *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gccuts.inc"
#include "geant321/gcjloc.inc"
#include "geant321/gconsp.inc"
#include "geant321/gcphys.inc"
#include "geant321/gcstak.inc"
#include "geant321/gctmed.inc"
#include "geant321/gctrak.inc"
#include "geant321/gcvolu.inc"
#include "geant321/gcvol1.inc"
#include "geant321/gcnum.inc"
#include "geant321/gcunit.inc"

*
* ** The following common is in GTMEDI. LSAMVL is set to true if
* ** we did not change volume yet
*
      COMMON/GCCHAN/LSAMVL
      LOGICAL LSAMVL
*
      DIMENSION R(3),D(3),U(3),QQ(3),vin(3)
      DIMENSION VBOU(3)
      LOGICAL LOLDTR
#if !defined(CERNLIB_SINGLE)
      PARAMETER (EPSMAC=1.E-6)
#endif
#if defined(CERNLIB_SINGLE)
      PARAMETER (EPSMAC=1.E-11)
#endif
      PARAMETER (MXPUSH=10)
      SAVE RIN1,EFFIC,DERI
C.
C.    ------------------------------------------------------------------
*
* *** Update local pointers if medium has changed
*
      LOLDTR=.FALSE.
      IF(IUPD.EQ.0)THEN
         IUPD  = 1
         JTCKOV = LQ(JTM-3)
         IF(JTCKOV.EQ.0) THEN
*
* *** This Cerenkov photon has crossed into a black medium.
* *** Just absorb it.
            IPROC = 101
            STEP  = 0.
            SLABS = 0.
            ISTOP = 2
            DESTEP = 0.
            GOTO 110
         ENDIF
         NPCKOV = Q(JTCKOV+1)
         JABSCO = LQ(JTCKOV-1)
         JEFFIC = LQ(JTCKOV-2)
         JINDEX = LQ(JTCKOV-3)
         JDERIV = LQ(JTCKOV-5)
         JPOLAR = LQ(JSTAK-1)
      ENDIF
      IF(SLENG.LE.0.) THEN
*
* *** Calculate GEKRAT for the particle
         IF(VECT(7).GE.Q(JTCKOV+NPCKOV+1)) THEN
            GEKRAT=1.
            IEKBIN=NPCKOV-1
         ELSEIF(VECT(7).LT.Q(JTCKOV+2)) THEN
*
* *** Particle below energy threshold ?  Short circuit
* *** This should never happen because the photons are generated
* *** only above threshold
*
*            GEKIN = 0.
*            GETOT = 0.
*            VECT(7)= 0.
*            ISTOP = 2
*            NMEC = 1
*            LMEC(1)= 30
*            GO TO 110
            GEKRAT=0.
            IEKBIN=1
         ELSE
            JMIN = 1
            JMAX = NPCKOV
   10       JMED = (JMIN+JMAX)/2
            IF(Q(JTCKOV+JMED+1).LT.VECT(7)) THEN
               JMIN = JMED
            ELSE
               JMAX = JMED
            ENDIF
            IF(JMAX-JMIN.GT.1) GO TO 10
            IEKBIN = JMIN
            GEKRAT = (VECT(7) - Q(JTCKOV+IEKBIN+1))/
     +      (Q(JTCKOV+IEKBIN+2)-Q(JTCKOV+IEKBIN+1))
         ENDIF
         GEKRT1=1.-GEKRAT
         RIN1=Q(JINDEX+IEKBIN)*GEKRT1+Q(JINDEX+IEKBIN+1)*GEKRAT
         DERI=Q(JDERIV+IEKBIN)*GEKRT1+Q(JDERIV+IEKBIN+1)*GEKRAT
         EFFIC=Q(JEFFIC+IEKBIN)*GEKRT1+Q(JEFFIC+IEKBIN+1)*GEKRAT
         STEPLA=Q(JABSCO+IEKBIN)*GEKRT1+Q(JABSCO+IEKBIN+1)*GEKRAT
      ENDIF
*
* *** Compute current step size
*
      IPROC  = 103
      STEP   = STEMAX
*     
*  **   Step limitation due to in flight absorbtion ?
*
      IF (ILABS.GT.0) THEN
         SLABS  = STEPLA*ZINTLA
         IF (SLABS.LT.STEP) THEN
            STEP  = SLABS
            IPROC = 101
         ENDIF
      ENDIF
*
      IF (STEP.LT.0.) STEP = 0.
*
*  **   Step limitation due to geometry ?
*
      STEPT=0.
      IF (STEP.GE.SAFETY) THEN
         CALL GTNEXT
         IF (IGNEXT.NE.0) THEN
*
* **    We are going to cross a boundary, so we need to simulate
* **    boundary effects and to know what is on the other side.
* **    For the moment save the current vector in the geometry tree.
*

            call gtmany(1)
*
* *** This is different from the other tracking routines.
* *** We get to the boundary and then we just jump over it
* *** So, linear transport till we are very near the boundary
*
#if !defined(CERNLIB_IBM)
            STEP = MAX(SNEXT-PREC,0.)
#endif
#if defined(CERNLIB_IBM)
            STEP = MAX(SNEXT-2.*PREC,0.)
#endif
C
C In case of Cherenkovs beeing near a corner, particle holding is
C prevented by a bigger step size. The value of VBOU is only used to
C calculate the correct surface normal (by GLISUR and GGPERP). The
C tracking of the Cherenkov is still using STEP
C
            IF(SNEXT.GT.0.) THEN
               DO 25 I=1,3
                  VBOU(I)=VECT(I)+SNEXT*VECT(I+3)
   25          CONTINUE
            END IF
C
            IF(STEP.GT.0.) THEN
               DO 30 I=1,3
                  VECT(I)=VECT(I)+STEP*VECT(I+3)
   30          CONTINUE
            ENDIF
            STEPT=STEP
#if !defined(CERNLIB_IBM)
            STEP  = SNEXT - STEP + PREC
#endif
#if defined(CERNLIB_IBM)
            STEP  = SNEXT - STEP + 2.*PREC
#endif
            IPROC = 0
            INWVOL= 2
            NMEC =  1
            LMEC(1)=1
         ENDIF
*
*        Update SAFETY in stack companions, if any
*        This may well not work.
         IF (IQ(JSTAK+3).NE.0) THEN
            DO 40 IST = IQ(JSTAK+3),IQ(JSTAK+1)
               JST = JSTAK +3 +(IST-1)*NWSTAK
               Q(JST+11) = SAFETY
   40       CONTINUE
            IQ(JSTAK+3) = 0
         ENDIF
*
      ELSE
         IQ(JSTAK+3) = 0
      ENDIF
*
* *** Linear transport
*
      VIN(1) = VECT(1)
      VIN(2) = VECT(2)
      VIN(3) = VECT(3)
      
      IF (INWVOL.EQ.2) THEN
         NBPUSH = 0
   50    DO 60 I = 1,3
            VECTMP  = VECT(I) +STEP*VECT(I+3)
            IF(VECTMP.EQ.VECT(I)) THEN
*
* *** Correct for machine precision
*
               IF(VECT(I+3).NE.0.) THEN
                  VECTMP = VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*
     +            EPSMAC
                  IF(NMEC.GT.0) THEN
                     IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
                  ENDIF
                  NMEC=NMEC+1
                  LMEC(NMEC)=104
#if defined(CERNLIB_DEBUG)
                  WRITE(CHMAIL, 10000)
                  CALL GMAIL(0,0)
                  WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
                  CALL GMAIL(0,0)
10000 FORMAT(' Boundary correction in GTCKOV: ',
     +       '    GEKIN      NUMED       STEP      SNEXT')
10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
#endif
               ENDIF
            ENDIF
            VOUT(I) = VECTMP
   60    CONTINUE
         NUMED1=NUMED
         CALL GTMEDI(VOUT,NUMED2)
         LOLDTR=.TRUE.
         IF(NUMED2.LE.0) THEN
            VECT(1) = VOUT(1)
            VECT(2) = VOUT(2)
            VECT(3) = VOUT(3)
            GO TO 110
         ENDIF
         ISH = 0
         IF(JVOLUM.GT.0) THEN
            JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
            ISH = Q(JVO+2)
         ENDIF    
         IF(LSAMVL.AND.ISH.NE.12) THEN
*
* *** In spite of our efforts we have not crossed the boundary
* *** we increase the step size and try again
*
            NBPUSH = NBPUSH + 1
            IF (NBPUSH.LE.MXPUSH) THEN
              STEP = STEP + NBPUSH*PREC
              GOTO 50
            ELSE
              INWVOL = 0
            ENDIF

         ENDIF
         IF(NUMED1.EQ.NUMED2) THEN
*
* *** If we are in the same medium, nothing needs to be done!
*
            VECT(1)=VOUT(1)
            VECT(2)=VOUT(2)
            VECT(3)=VOUT(3)
            IPROC=0
         ELSE
            JTM2 = LQ(JTMED-NUMED2)
            IF(IQ(JTM2-2).GE.3) THEN
               JTCKV2 = LQ(JTM2-3)
            ELSE
               JTCKV2 = 0
            ENDIF
            IF(JTCKV2.GT.0) THEN
               NPCKV2 = Q(JTCKV2+1)
               JABSC2 = LQ(JTCKV2-1)
               JEFFI2 = LQ(JTCKV2-2)
               JINDX2 = LQ(JTCKV2-3)
               IF(VECT(7).GE.Q(JTCKV2+NPCKV2+1)) THEN
                  GEKRT2=1.
                  IEKBI2=NPCKV2-1
               ELSEIF(VECT(7).LT.Q(JTCKV2+2)) THEN
                  GEKRT2=0.
                  IEKBI2=1
               ELSE
                  JMIN = 1
                  JMAX = NPCKV2
   64             JMED = (JMIN+JMAX)/2
                  IF(Q(JTCKV2+JMED+1).LT.VECT(7)) THEN
                     JMIN = JMED
                  ELSE
                     JMAX = JMED
                  ENDIF
                  IF(JMAX-JMIN.GT.1) GO TO 64
                  IEKBI2 = JMIN
                  GEKRT2 = (VECT(7) - Q(JTCKV2+IEKBI2+1))/
     +            (Q(JTCKV2+IEKBI2+2)-Q(JTCKV2+IEKBI2+1))
               ENDIF
               GEKRT1=1.-GEKRT2
               ABSCO2=Q(JABSC2+IEKBI2)*GEKRT1+Q(JABSC2+IEKBI2+1)*GEKRT2
               EFFIC2=Q(JEFFI2+IEKBI2)*GEKRT1+Q(JEFFI2+IEKBI2+1)*GEKRT2

               IF(JINDX2.GT.0) THEN
                 RIN2=Q(JINDX2+IEKBI2)*GEKRT1+Q(JINDX2+IEKBI2+1)*GEKRT2
               ELSE
                 RIN2=0.
               ENDIF
               IF(ABS(RIN2-RIN1).LE.PREC) THEN
*
* *** The two media have the same refraction coefficient, i.e. nothing
* *** has to be done
*
                  VECT(1)=VOUT(1)
                  VECT(2)=VOUT(2)
                  VECT(3)=VOUT(3)
                  STEPLA = ABSCO2
                  EFFIC  = EFFIC2
                  IPROC=0
               ELSE
                  IPROC = 102
               ENDIF
            ELSE
               ISTOP=2
               DESTEP=0
               NMEC = NMEC+1
               LMEC(NMEC)=30
               VECT(1)=VOUT(1)
               VECT(2)=VOUT(2)
               VECT(3)=VOUT(3)
               GOTO 71
            ENDIF
         ENDIF
      ELSE
         DO 70 I = 1,3
            VECT(I)  = VECT(I) +STEP*VECT(I+3)
   70    CONTINUE
      ENDIF
*
   71 STEP = STEPT + STEP
      SLENG = SLENG +STEP
*
* *** Update time of flight 
* *** using the group velocity (AM)
*
      TOFG = TOFG +STEP/CLIGHT*RIN1*(1.+DERI) 
*
* *** Update interaction probabilities
*
      IF (ILABS.GT.0)    ZINTLA = ZINTLA -STEP/STEPLA
*
      IF (IPROC.EQ.0) GO TO 110
      NMEC = NMEC+1
      LMEC(NMEC) = IPROC
*
*  ** Absorbtion in flight ?
*
      IF (IPROC.EQ.101) THEN
         ISTOP=2

         CALL GRNDM(RNDM,1)
         IF(RNDM.LT.EFFIC) THEN
*
* ***  Destep =/= 0 means that the photon has been detected
*

            DESTEP=VECT(7)
         ELSE
            DESTEP=0.
         ENDIF
*
*  ** Boundary action?
*
      ELSE IF (IPROC.EQ.102) THEN
         IF(JINDX2.EQ.0) THEN
*
* *** Case dielectric -> metal
*
            CALL GRNDM(RNDM,1)
            IF(RNDM.LT.ABSCO2) THEN
*
* *** Photon is absorbed in the next medium
*
               NMEC=NMEC+1
               LMEC(NMEC)=101
               ISTOP = 2
               CALL GRNDM(RNDM,1)
               IF(RNDM.LT.EFFIC2) THEN
*
* ***  Destep =/= 0 means that the photon has been detected
*
                  DESTEP=VECT(7)
               ELSE
                  DESTEP = 0.
               END IF
               VECT(1) = VOUT(1)
               VECT(2) = VOUT(2)
               VECT(3) = VOUT(3)
               GOTO 110
            ELSE
*
* *** Photon is reflected (no polarization effects)
*
               CALL G3LISUR(VECT,VBOU,NUMED1,NUMED2,U,PDOTU,IERR)
               IF (IERR.NE.0) THEN
                  WRITE(CHMAIL,10200) IERR
                  CALL GMAIL(0,0)
                  GO TO 110
               END IF
*
* ** Restore old volume tree, the photon does not cross the boundary
*
*              CALL GLVOLU(NLEVL1,NAMES1,NUMBR1,IERR)
               CALL GTMEDI(VIN,N)
               LOLDTR=.FALSE.
*
*
* *** Reflect, but make sure we are in the old volume before
*
               NBPUSH = 0
   80          NBPUSH = NBPUSH+1
               IF(NBPUSH.GT.MXPUSH) THEN
                  WRITE(CHMAIL,10300) NTMULT,NSTEP
                  CALL GMAIL(0,0)
                  ISTOP=1
                  GOTO 110
               ELSE
                  CALL GINVOL(VECT,ISAME)
                  IF(ISAME.EQ.0) THEN
                     PRECN = NBPUSH*PREC
                     VECT(1) = VECT(1) - PRECN*VECT(4)
                     VECT(2) = VECT(2) - PRECN*VECT(5)
                     VECT(3) = VECT(3) - PRECN*VECT(6)
                     GO TO 80
                  ENDIF
               ENDIF
*
               NMEC=NMEC+1
               LMEC(NMEC)=106
               EDOTU = POLAR(1)*U(1) + POLAR(2)*U(2) + POLAR(3)*U(3)
               VECT(4) = +VECT(4) - 2*PDOTU*U(1)
               VECT(5) = +VECT(5) - 2*PDOTU*U(2)
               VECT(6) = +VECT(6) - 2*PDOTU*U(3)
               POLAR(1) = -POLAR(1) + 2*EDOTU*U(1)
               POLAR(2) = -POLAR(2) + 2*EDOTU*U(2)
               POLAR(3) = -POLAR(3) + 2*EDOTU*U(3)
               INWVOL  = 0
            ENDIF
         ELSE
*
*  case dielectric-dielectric:
*
            CALL G3LISUR(VECT,VBOU,NUMED1,NUMED2,U,PDOTU,IERR)
            IF (IERR.NE.0) THEN
               WRITE(CHMAIL,10200) IERR
               CALL GMAIL(0,0)
               GO TO 110
            END IF
            EDOTU = POLAR(1)*U(1) + POLAR(2)*U(2) + POLAR(3)*U(3)
            COST1 = -PDOTU
            IF (ABS(COST1).LT.1.) THEN
               SINT1 = SQRT((1-COST1)*(1+COST1))
               SINT2 = SINT1*RIN1/RIN2
            ELSE
               SINT1 = 0.0
               SINT2 = 0.0
            END IF
            IF (SINT2.GE.1) THEN
*
* ***  Simulate total internal reflection
*
               NMEC=NMEC+1
               LMEC(NMEC)=106
*
* ** Restore old volume tree, the photon does not cross the boundary
*
*              CALL GLVOLU(NLEVL1,NAMES1,NUMBR1,IERR)
               CALL GTMEDI(VIN,N)
               LOLDTR=.FALSE.
*
* *** Reflect, but make sure we are in the old volume before
*
               NBPUSH = 0
   90          NBPUSH = NBPUSH+1
               IF(NBPUSH.GT.MXPUSH) THEN
                  WRITE(CHMAIL,10300) NTMULT,NSTEP
                  CALL GMAIL(0,0)
                  ISTOP=1
                  GOTO 110
               ELSE
                  CALL GINVOL(VECT,ISAME)
                  IF(ISAME.EQ.0) THEN
                     PRECN = NBPUSH*PREC
                     VECT(1) = VECT(1) - PRECN*VECT(4)
                     VECT(2) = VECT(2) - PRECN*VECT(5)
                     VECT(3) = VECT(3) - PRECN*VECT(6)
                     GO TO 90
                  ENDIF
               ENDIF
*
               VECT(4) = +VECT(4) - 2*PDOTU*U(1)
               VECT(5) = +VECT(5) - 2*PDOTU*U(2)
               VECT(6) = +VECT(6) - 2*PDOTU*U(3)
               POLAR(1) = -POLAR(1) + 2*EDOTU*U(1)
               POLAR(2) = -POLAR(2) + 2*EDOTU*U(2)
               POLAR(3) = -POLAR(3) + 2*EDOTU*U(3)
               INWVOL = 0
            ELSE
*
* ***  Calculate amplitude for transmission (Q = P x U)
*
               COST2 = SIGN(1.0,COST1)*SQRT((1-SINT2)*(1+SINT2))
               IF (SINT1.GT.1.E-4) THEN
                  QQ(1) = VECT(5)*U(3) - VECT(6)*U(2)
                  QQ(2) = VECT(6)*U(1) - VECT(4)*U(3)
                  QQ(3) = VECT(4)*U(2) - VECT(5)*U(1)
                  EPERP1 = (POLAR(1)*QQ(1) + POLAR(2)*QQ(2) + POLAR(3)*
     +            QQ(3))
                  ENORM  = SQRT(EPERP1**2+EDOTU**2)
                  EPERP1 = EPERP1/ENORM
                  EPARL1 = EDOTU/ENORM
               ELSE
                  QQ(1) = POLAR(1)
                  QQ(2) = POLAR(2)
                  QQ(3) = POLAR(3)
*
*     Here we follow Jackson's conventions and we set the parallel
*     component = 1 in case of a ray perpendicular to the surface
                  EPERP1 = 0.
                  EPARL1 = 1.
               END IF
               IF(COST1.NE.0.) THEN
                  S1 = RIN1*COST1
                  EPERP2 = 2*S1*EPERP1/(RIN1*COST1+RIN2*COST2)
                  EPARL2 = 2*S1*EPARL1/(RIN2*COST1+RIN1*COST2)
                  E2 = EPERP2**2 + EPARL2**2
                  S2 = RIN2*COST2*E2
                  TCOEF = S2/S1
               ELSE
                  TCOEF = 0.
               ENDIF
               CALL GRNDM(RNDM,1)
               IF (RNDM.GT.TCOEF) THEN
*
* *** Simulate reflection
*
                  NMEC=NMEC+1
                  LMEC(NMEC)=106
*
* *** Restore old volume tree, the photon does not cross the boundary
*
*                 CALL GLVOLU(NLEVL1,NAMES1,NUMBR1,IERR)
                  CALL GTMEDI(VIN,N)
                  LOLDTR=.FALSE.
*
* *** Reflect, but make sure we are in the old volume before
*
                  NBPUSH = 0
  100             NBPUSH = NBPUSH+1
                  IF(NBPUSH.GT.MXPUSH) THEN
                     WRITE(CHMAIL,10300) NTMULT,NSTEP
                     CALL GMAIL(0,0)
                     ISTOP=1
                     GOTO 110
                  ELSE
                     CALL GINVOL(VECT,ISAME)
                     IF(ISAME.EQ.0) THEN
                        PRECN = NBPUSH*PREC
                        VECT(1) = VECT(1) - PRECN*VECT(4)
                        VECT(2) = VECT(2) - PRECN*VECT(5)
                        VECT(3) = VECT(3) - PRECN*VECT(6)
                        GO TO 100
                     ENDIF
                  ENDIF
*
                  VECT(4) = VECT(4) - 2*PDOTU*U(1)
                  VECT(5) = VECT(5) - 2*PDOTU*U(2)
                  VECT(6) = VECT(6) - 2*PDOTU*U(3)
                  IF(SINT1.GT.1E-4) THEN
                     EPARL2 = RIN2*EPARL2/RIN1-EPARL1
                     EPERP2 = EPERP2-EPERP1
                     E2 = EPERP2**2 + EPARL2**2
                     R(1) = U(1) + PDOTU*VECT(4)
                     R(2) = U(2) + PDOTU*VECT(5)
                     R(3) = U(3) + PDOTU*VECT(6)
                     EABS = SQRT(E2)*SINT1
                     CPARL = EPARL2/EABS
                     CPERP = EPERP2/EABS
                     POLAR(1) = CPARL*R(1) - CPERP*QQ(1)
                     POLAR(2) = CPARL*R(2) - CPERP*QQ(2)
                     POLAR(3) = CPARL*R(3) - CPERP*QQ(3)
                  ELSEIF(RIN2.GT.RIN1) THEN
*
* *** Case of ray perpendicular to the surface. No change or
* *** an inversion of phase.
                     POLAR(1) = -POLAR(1)
                     POLAR(2) = -POLAR(2)
                     POLAR(3) = -POLAR(3)
                  ENDIF
                  INWVOL = 0
               ELSE
*
* *** Simulate transmission/refraction
*
                  NMEC=NMEC+1
                  LMEC(NMEC)=107
                  VECT(1) = VOUT(1)
                  VECT(2) = VOUT(2)
                  VECT(3) = VOUT(3)
                  IF(SINT1.GT.1E-4) THEN
                     ALPHA = COST1-COST2*(RIN2/RIN1)
                     D(1) = VECT(4) + ALPHA*U(1)
                     D(2) = VECT(5) + ALPHA*U(2)
                     D(3) = VECT(6) + ALPHA*U(3)
                     DABS = SQRT(D(1)**2+D(2)**2+D(3)**2)
                     VECT(4) = D(1)/DABS
                     VECT(5) = D(2)/DABS
                     VECT(6) = D(3)/DABS
                     PDOTU = -COST2
                     R(1) = U(1) - PDOTU*VECT(4)
                     R(2) = U(2) - PDOTU*VECT(5)
                     R(3) = U(3) - PDOTU*VECT(6)
                     EABS = SQRT(E2)
                     CPARL = EPARL2/(EABS*SINT2)
                     CPERP = EPERP2/(EABS*SINT1)
                     POLAR(1) = CPARL*R(1) + CPERP*QQ(1)
                     POLAR(2) = CPARL*R(2) + CPERP*QQ(2)
                     POLAR(3) = CPARL*R(3) + CPERP*QQ(3)
                  ENDIF
                  GEKRAT = GEKRT2
                  IEKBIN = IEKBI2
                  STEPLA = ABSCO2
                  EFFIC = EFFIC2
                  RIN1 = RIN2
               END IF
            END IF
         END IF
      ENDIF
*                                                             END GTCKOV
  110 IF(LOLDTR) CALL GLVOLU(NLEVL1,NAMES1,NUMBR1,IERR)
10200  FORMAT(' **** GTCKOV: error from GLISUR = ',I10)
10300  FORMAT(' **** GTCKOV: unable to reflect at NTMULT = ',
     +        I8,' step No. ',I8,' photon abandoned!')
      END
